Imaging local diffusion in microstructures using NV-based pulsed field gradient NMR

Understanding diffusion in microstructures plays a crucial role in many scientific fields, including neuroscience, medicine, or energy research. While magnetic resonance (MR) methods are the gold standard for diffusion measurements, spatial encoding in MR imaging has limitations. Here, we introduce nitrogen-vacancy (NV) center–based nuclear MR (NMR) spectroscopy as a powerful tool to probe diffusion within microscopic sample volumes. We have developed an experimental scheme that combines pulsed gradient spin echo (PGSE) with optically detected NV-NMR spectroscopy, allowing local quantification of molecular diffusion and flow. We demonstrate correlated optical imaging with spatially resolved PGSE NV-NMR experiments probing anisotropic water diffusion within an individual model microstructure. Our optically detected PGSE NV-NMR technique opens up prospects for extending the current capabilities of investigating diffusion processes with the future potential of probing single cells, tissue microstructures, or ion mobility in thin film materials for battery applications.


Short summary of the working principle of the NV center
The NV center is a point defect in the diamond crystal lattice (Fig. S1b), where one of the carbon atoms is replaced by a nitrogen atom and an adjacent lattice site is left vacant, resulting in a C 3v symmetry. The NV center can exists in different charge states, but only the negatively charged NV − state is applicable for quantum sensing. Therefore, we will refer to it as "NV center". The NV center's electronic ground and excited states are triplet states, located within the diamond band-gap (64). The corresponding m s = 0 and m s = ± 1 in the electronic ground state are separated by a zero-field splitting of D = 2.87 GHz, whereas the m s = ± 1 degeneracy can be lifted according to the Zeeman splitting ∆E = 2γB, where γ ≈ 28 GHz/T is the electron's gyromagnetic ratio. The NV-center can be excited with green light (green arrows, Fig. S1a) and will decay back into its respective electronic ground state by emitting red fluorescence (red arrows, Fig. S1a). The m s = ± 1 spin states have an additional, non-radiative decay channel via singlet states, which in turn can decay into the m s = ± 1 as well as m s = 0 states. This leads to i) the possibility of initializing the NV center in the m s = 0 state by optical pumping, and ii) to a spin-state dependent fluorescence, which allows to read out the NV-center's quantum state optically. Lastly, the magnetic sublevels can be coherently addressed with microwave pulses (turquoise arrows, Fig. S1a).

NV-AC magnetometry using CASR
For most AC-sensing applications (such as the detection of NMR signals) the NV-center is dynamically decoupled from the environment by a train of π-pulses with a spacing 2τ . In the simplest case, this is the Carr-Purcell-Meiboom-Gill (CPMG) sequence (61,65). We use one of its variations called universally robust dynamic decoupling (UDD) (58), as it is more robust against pulse errors. These pulse sequences act as a narrow-band frequency filter around f sense = 1 4τ (Fig. S2): Magnetic field frequencies not matching this condition will be suppressed by the spin echoes of the dynamic decoupling pulse sequence. Frequencies matching this condition will lead to a phase accumulation on the NV-center's spin state (related to electron spin echo envelope modulation "ESEEM" in conventional spin resonance experiments), which can be read out optically and is used herein for the detection of NMR signals. In our experiments we use a modified version called coherently averaged synchronized readout (CASR) (24). It relies on the fact that the signal detected by dynamical decoupling sequences is sensitive to the relative phase of the ac signal (NMR signal) and the pulse sequence (i.e., the train of π pulses).
If the target ac signal is asymmetric around the π-pulse, then the phase accumulation is maximal. This is due to the simultaneous "change in sign" of both the oscillating magnetic field and the NV-center's spin. If in turn the ac signal is symmetric around the π-pulses, no phase will be acquired and the quantum state of the NV is not influenced. If the period of the sensing frequency is an integer multiple of each measurement repetition, the phase relation will be constant between consecutive readouts. However, if the frequency is not matched as depicted in Fig. S2, we observe aliasing and the fluorescence readout will oscillate at the difference frequency. The advantage of this pulse sequence is, that NMR signals can be detected with high spectral resolution, since the only limitation is the lifetime of the NMR-active nuclei under investigation (assuming insignificant back-action between the NV and samples spins). A more detailed description of this process can be found in the literature (24,56,63,66). Figure S2: Detecting NMR signals using the CASR pulse sequence. A sketch of CASR's working principle. The spacing of the individual measurement blocks (dynamic decoupling sequence and the optical readout pulse) is a multiple of a frequency close to the target ACsignal (blue). The difference in the initial phase of the AC-signal relative to the start of each individual block (green rectangles), leads to an oscillation of the phase accumulation of the NVcenter's spin state over consecutive measurements. This, in turn, causes an oscillation of the NV fluorescence intensity in the recorded data (red). This is illustrated by the two insets in the bottom right, the first having a relative phase leading to minimal signal amplitude, the second leading to maximum signal amplitude.   Figure S5: Photos of the experimental setup. a) Superconducting magnet mounted on a tiltable and rotating platform for magnetic field alignment. The probe head is mounted on a aluminum rod. b) Photo of the probe head. The two lenses (II) are used for optical imaging, while the excitation laser is focused on the diamond chip using lens I. c) The 3D-printed sample holder. The excitation laser is reflected onto the diamond chip using a mirror at I, while the diamond is imaged using mirror II. The magnetic field gradient coils are stacked into the holder (the dark green PCBs). d) A zoom-in on the diamond chip. The copper-coils in Helmholtz-geometry are used for driving the sample nuclear spins. A small RF coil (right) is used for calibrating the experiment. A shortened coaxial cable used for driving the NV centers is mounted on the microfluidic chip (center).

PGSE in a microfluidic channel
In the following section, we outline the basic theory of restricted diffusion and discuss the modeling of our diffusion experiments. The theory behind sample-diffusion is well understood (6) including the case of restricted diffusion on the microscale (67-69).
A Analytical treatment and numerical simulation of restricted diffusion on short time scales in one dimension, e.g. between two parallel plates, can be calculated according to: where L is the distance between the two plates x and (x 0 − x (t)) 2 is the mean squared particle displacement. These results are identical to the diffusion between two parallel walls of a cuboid microfluidic channel, and represent the eigenvalues of the diffusion tensor. In the case of the microfluidic channel this would be the basis, in which each of the basis vectors is orthogonal to one set of the channel-walls. The measurement basis in our experiments follows the NMRconvention (B 0 ||ẑ), which is equivalent to a rotation of the previous coordinate system around x by φ ≈ 35.26 • , withx =x,ŷ →ŷ andẑ →ẑ. This leads to the effective diffusion constants of: .

(S2)
The same results can be obtained by numerical simulations using a random walk in a sample volume of the size of our microfluidic channel. At each time step ∆t, each sample point is moved by a randomly Gaussian distributed amount. If it would leave the sample volume, it is reflected back inwards. From this we calculated the mean-squared displacement and in turn the effective diffusion constant along each direction, according to: A comparison of the analytical and numerically determined effective diffusion constant is shown in Fig. S7.

B Influence of the NV center's sensitivity map on diffusion measurements
In this paragraph, we discuss the possible influence of the strongly geometry-dependent NV-NMR detection on the PGSE-experiments (30). The field of a magnetic dipole, e.g. a nuclear spin in our sample, is described by: where µ 0 is the magnetic permeability,r is the direction from the sample spin to the NV center, and r = | r | is the distance between the two. m(t) is the magnetic moment of the precessing nuclear spin. NV-NMR typically measures the parallel magnetic field component, via the Zeeman shift of the NV center's |±1 levels. Projecting B spin ( r, t) onto B 0 , we obtain: Since the magnetization of the precessing spins are orthogonal toB 0 , after a π/2-pulse, ( m(t) ·B 0 ) = 0, leading to a magnetic field contribution B spin ( r, t) of a sample-spin at position r: A visualization of this sensitivity map can be found in Fig. S9 a) and b). In particular the two terms 1/r 3 and (r ·B 0 ) can have a strong impact on the directional sensitivity of our NV centers. Since the time independent contribution to the signal amplitude is of importance, we integrate over one period of sample spin oscillation T and arrive at a relative position dependence c of each individual spin, according to: where C is chosen such that c( r) dV sample = 1 .
We note that, for our purposes, this integral is convergent as a function of sample volume (30), therefore the bulk of the signal stems from a local area close to the diamond surface, whose size is of the same order of magnitude as the depth of the NV-center. For highly anisotropic samples, an additional contribution to the signal can appear, since sample areas with positive and negative signal contributions at the NV-location might not cancel out. Nevertheless, even in a worst-case scenario, the total signal from dipoles with a distance from the NV center much bigger than the NV center's depth below the diamond surface will only increase logarithmically as a function of sample volume (30). If the NV-layer is chosen such that the layer thickness is not several orders of magnitude larger than the smallest sample dimension, this additional signal can probably be considered to be insignificant.
Finally, we can define a position-weighted diffusion coefficient in each direction according to: The diffusion coefficient weighted by the coupling to an ensemble of NVs, distributed over a cylindrical volume of 45 µm diameter and 50 µm depth was simulated. The resulting difference between the weighted diffusion coefficients with respect to their analytical solutions from equation S2 is plotted in Fig. S10. The difference between the weighted simulation and analytical results stays below 1 % of the analytical value for ∆ = 0.5 s in theŷ direction and below 1 % for bothx andẑ for more than ∆ = 0.5 s. In the interpretation of our data, this difference is therefore treated as negligible.
We also note that, although a significant fraction of the molecules may diffuse out of the detection volume and be replaced by other molecules after the second gradient pulse, this will not affect the recorded signal. We assume that both the B0 field and the RF field, driving the nuclear sample spins, are homogeneous over an area much larger than i) the size of the laser spot, which defines the detection volume, and ii) the mean distance diffused during the measurement time. For example, the sample spins probed in the different locations in Fig.   5 in the main manuscript all have indistinguishable Rabi frequencies, despite their locations are several hundred micrometers apart. Assuming the diffusion constant and environment, i.e.
local anisotropicity, of all spins entering and leaving the probed volume is the same, they should underlie the same Gaussian process and lead to the same signal attenuation, resulting in no net effect on the recorded diffusion constant. Therefore, the so-called "diffusion blurring" during the diffusion weighting with magnetic field gradients will not contribute to the reduction of the spatial resolution. Contrary to that, diffusion processes during the acquisition of the spin-echo from the sample will contribute with the conventional square-root dependence of the induced blurring on the acquisition time, which mandates using short data acquisition windows if high spatial resolution is desired. Figure S10: Impact of NV center sensitivity map on diffusion measurements. The difference in percent between the diffusion tensors weighted diagonal elements with respect to the unweighted analytical value as a function of gradient pulse spacing ∆.

Estimating the impact of a heterogeneous flow profile
The laminar flow in our microfluidic channel has a quadratic velocity distribution across the cross-section of the channel. In the case of a homogeneous sensitivity over the whole channel, the phase-shift measured will depend on the average velocity in the channel, which is set as constant at the flow pump and therefore independent of any velocity distribution. This means the only effect of a quadratic flow profile should be signal attenuation. However, the assumption of a spatially homogeneous sensitivity does not apply for NV-NMR (see section 4).
We used numerical simulations to analyze the effect of the interplay between the quadratic flow profile and the sensitivity map of the NV center (Fig. S11). The maximum velocity of the quadratic flow profile v max relates to the set or average velocity v mean as: v max = 3/2 v mean .
Since the phase shift, δφ induced in the PGSE experiment is linear with the velocity, the flow profile was included in the simulation using a quadratic phase shift over the channel's cross section with a maximum phase shift δφ max = 3/2 δφ mean . For 201 linearly spaced points between δφ mean = 0 and δφ mean = 2π the corresponding NV-NMR signal phase δφ simulated was simulated for ∆ = 0.3 s. As shown in Fig. S11, the simulated signal phase shift δφ simulated initially deviates slightly from the expected linear dependence on the mean sample phase shift δφ mean , though for δφ mean 1.5 π the deviation becomes significant. In the experiments performed in the main manuscript, the expected mean sample phase shift does not exceed δφ mean ≈ 1.27 π. The discrepancy between mean sample phase shift and signal phase shift is corrected in the data displayed in the main manuscript by fitting a polynomial of order 4 to δφ mean (δφ simulated ). This fit maps the recorded measurement value to the expected mean velocity according to our numerical simulation. Figure S11: Simulated phase shift for velocimetry experiments. The extracted NV-NMR signal phase shift is plotted as a function of the set, mean phase shift (light blue). A linear dependency of δφ mean α δφ set is depicted with a dashed, dark blue line. Dashed, grey guiding lines were added for both δφ set = π and δφ mean = π.

Outlook and feasibility of high spatial resolution PGSE NV-NMR
Increasing the spatial resolution of PGSE NV-NMR by optical means will reduce the signal to noise (SNR) ratio due to a decrease in the number of interrogated NV centers. For instance, increasing the spatial resolution by a factor of 2, will lead to a reduction of ∼ 2 3 in the number of NV centers contributing to the detected signal (and the corresponding increase in averag-ing time). We assume the following parameters for the lowest requirement for a PGSE NV-NMR experiment: six different gradient strength steps, averaged until a maximal SNR of 10 is reached (without an applied magnetic field gradient), with a duration of t avr = 3 s for a single average. In Fig. S12

Laser and microwave induced sample heating
Since the ADC is temperature dependent (43), measuring the sample's temperature is important -in particular for NV-based experiments, where microwave or laser-induced heating can occur.
However, measuring temperature in microfluidic channels is challenging due to small sample volumes. We utilized the NV center's ability to measure temperature, as the energy splitting of its two-level system is temperature-dependent (70). Under our experimental conditions, we are only able the probe a temperature change as a function of microwave and laser power.
Moreover, the measured values give information on the temperature of the diamond -not of the sample (both are in thermal contact). We also note that the actual temperature may vary, depending on the experimental situation (e.g. position of the microwave loop, distance of the measurement to the microwave loop etc.). Therefore, these values should be considered as a rough estimate. Figure S13: Measurement of laser and MW induced sample heating. a) The pulse sequence used for measuring the sample heating. A chain of ESR-experiments is used to detect the relative temperature of the diamond. Before the experiment a MW or laser pulse of length 0.3 s is used to heat the sample. Each individual ODMR spectrum is fitted, and the temperature during this spectrum is deduced from its midpoint. b) The relative temperature drift during the PGSE sequence as a function of time, as measured by the NV centers for various MW frequencies and powers tested.
In our experiment, we repeated the PGSE sequence, substituting the universally robust dynamical decoupling sequence with ESR sequences to determine the resonance frequency of the NV-ensemble. The MW frequency is swept between f 0 ± 5 MHz, where f 0 is the resonance frequency of the NV centers and each pair of dips in NV-fluorescence corresponds to a single ODMR spectrum of the hyperfine splitting of the NV-centre's |0 → |−1 transition due to the 15 N-nucleus. Each individual spectrum is fit to a double Lorentzian model and f 0 is extracted.
For each dataset, the mean of the last five datapoints is subtracted, to cancel out the effects of small temporal shifts in the magnetic field. Due to technical constraints we were only able to measure the relative change in temperature and not the absolute temperature.
To check microwave induced heating, we vary the amplitude I and frequency f of a 0.3 s long MW-pulse before our experiment, corresponding to the ODNP pulse in the NV-NMRexperiments. Since this microwave pulse has by far the longest duration and highest duty cycle, it is expected to produce the most significant heating effects. Results of these measurements can be seen in Fig. S13. For the MW-frequencies used in our experiments, f ≈ 2.02 GHz and f ≈ 4.87 GHz, we observe no significant heating effects. Only at higher frequencies, e.g. at f ≈ 10.0 GHz, we see a significant heating effect of the NV-diamond.
To check for laser induced sample heating, the MW-pulse is replaced by a laser pulse. Here we observe small heating of ca. 0.5 K for a laser duty cycle of 100 % and no change when reducing the lasers duty cycle to 0%.
Since in our PGSE NV-NMR experiments much lower duty-cycles and/or amplitudes of both the laser and the MW are used, we exclude any excessive heating effects on our sample. Nevertheless, we anticipate minor heating that can influence the observed ADCs, as discussed in the main text.
Similar temperature calibration measurements as the ones done in this section could be performed by recording 1 H NMR spectra of ethylene glycol. However, only temperature drifts on a time scale longer than an NMR experiment could be recorded and fast temperature fluctuations induced by laser or microwave pulses would be invisible.

Individual ADC measurements
The ADC-tensors as calculated following Kingsley et al. (37).  10 Example data set of a measured NV-NMR PGSE experiment An example data set of spin-echos are shown in Fig. S15, both in the time as well as in the frequency domain. Figure S15: Example measured dataset of an NV-NMR PGSE experiment. A set of acquired spin-echos is displayed, both in the time-domain (left) as well as in the frequency-domain (right). This data set was measured using theẑ-direction gradient coils at the location marked as "location 3" in the main manuscript. The current through the gradient coils was swept from -5 A (bottom) to 5 A (top). The increase and decrease of the signal amplitude is clearly visible.